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We develop a method to infer log-normal random fields from measurement data affected by 
Gaussian noise. The log-normal model is well suited to describe strictly positive signals with fluctu- 
ations whose amplitude varies over several orders of magnitude. We use the formalism of minimum 
Gibbs free energy to derive an algorithm that uses the signal's correlation structure to regularize 
the reconstruction. The correlation structure, described by the signal's power spectrum, is thereby 
reconstructed from the same data set. We further introduce a prior for the power spectrum that 
enforces spectral smoothness. The appropriateness of this prior in different scenarios is discussed 
and its effects on the reconstruction's results are demonstrated. We validate the performance of our 
reconstruction algorithm in a series of one- and two-dimensional test cases with varying degrees of 
non-linearity and different noise levels. 
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I. INTRODUCTION 

Reconstructing continuous fields from a finite and 
noisy data set is a problem encountered often in all 
branches of physics and the geo-sciences. In this paper 
we develop a method to reconstruct a log-normal field, 
i.e. a field whose logarithm can be modeled as a Gaussian 
random field, defined on an arbitrary manifold. We show 
simple one-dimensional examples as well as planar and 
spherical two-dimensional cases. 

The log-normal model is well suited to describe many 
physical fields. Its main features are that the field values 
are guaranteed to be positive, that they may vary across 
many orders of magnitude, and that spatial correlations 
may exist. One prominent example from the discipline 
of astrophysics that exhibits these three features and is 
therefore often modeled as a log-normal field is the mat- 
ter density in today's universe [see e.g. Ej- The simplest 
way to observationally estimate the matter density is to 
count the galaxies per volume. The relationship between 
the matter density and the galaxy number counts can 
be approximated as a Poisson process [e.g. HH [Mj , thus 
making the data likelihood Poissonian as well with the 
expected number of galaxies per volume element given 
as a function of the underlying log-normal density field. 
In the statistical literature, such a combination of a log- 
normal field with Poissonian statistics is known as a log 
Gaussian Cox process [T7]. This model has been applied 
in fields as diverse as finance [T], agriculture [3], atmo- 
spheric studies [2^ , and epidemiology [21 14] . 

The main motivation for our work, however, comes 
from astronomical observations of radiative processes. 
The intensity of radiation coming from different direc- 
tions in such observations can vary over many orders of 
magnitude while being spatially correlated. Thus ap- 
plying a log-normal model for its statistical description 
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seems natural. While photon counts bring with them a 
Poissonian likelihood, the observational uncertainty can 
be approximatively regarded to be Gaussian in the limit 
of large photon numbers. Therefore, we restrict ourselves 
to cases in which the measurement noise can be assumed 
to be additive and Gaussian. 

Apart from the noise contribution, we will assume a 
deterministic linear relationship between the log-normal 
field and the data. This model is general enough to ac- 
comodate a large variety of observational settings, such 
as targeted point-wise observations across limited areas, 
convolutions of the log-normal field, interferometric ob- 
servations leading to a Fourier transformed version of the 
log-normal field, or any combination of these effects. 

The inclusion of spatial correlations in the field model 
is necessary for an accurate statistical description and 
will aid in the reconstruction. If the spatial correlation 
structure is known, knowledge of the field values at some 
locations can be used to extrapolate the field into re- 
gions where the field value has not been measured or 
the observations have a higher uncertainty. However, in 
general it is not known a priori how strongly the field is 
correlated across a given distance. In the case of a sta- 
tistically homogeneous field, the correlation structure is 
described by the field's power spectrum. So if one wants 
to make use of the spatial correlations during the recon- 
struction, one needs to simultaneously infer the power 
spectrum. Several techniques have been developed to 
solve this problem for Gaussian fields [e.g [131 [26]. One 
such technique was derived within the framework of in- 
formation field theory \7\ [T^j in [B] and later rederived in 
[8], where the formalism of minimum Gibbs free energy 
was introduced and employed to tackle this problem. In 
the same paper, the problem of inferring a log-normal 
field with unknown power spectrum in the presence of 
Poissonian noise was briefly discussed. 

Here, we use the formalism of minimum Gibbs free 
energy to derive filter equations for the estimation of a 
log-normal field from data including a Gaussian noise 
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component. The filter we present results in an estimate 
of the field, its power spectrum, and the respective un- 
certainties. 

A problem that often arises when reconstructing power 
spectra is that the reconstructed power spectra do not 
turn out smooth, especially on scales with little redun- 
dancy, despite the true underlying spectra being smooth. 
One reason why this is not desirable is that sharp kinks 
and drop-offs in the power-spectrum can lead to ring- 
ing effects in the reconstructed field estimate. Another 
reason is that many of the fields occuring in nature are 
actually expected to have smoothly varying power over 
different scales since neighboring scales interact. One 
prominent example is the power spectrum of a viscous 
turbulent fluid [15 . Smoothness of the power spectrum is 
often enforced by an ad-hoc smoothing step in the power 
spectrum reconstruction [l8l - [20] . Here, we follow an idea 
presented in [5" and show how to enforce spectral smooth- 
ness by introducing an appropriate prior for the power 
spectrum that punishes non-smooth spectra. We demon- 
strate the feasibility of this approach using one specific 
example for such a smoothness prior both for the recon- 
struction of Gaussian and log-normal fields and discuss 
the range of applicability of the chosen spectral smooth- 
ness prior. Our approach also allows for the estimation 
of a complete uncertainty matrix for the estimated power 
spectrum. 

We first develop the formalism of the spectral smooth- 
ness prior for the case of Gaussian signal fields in Sec.|TTj 
where we show how to derive the filter formulas of [6j 
and [8] in a different way that easily accomodates an 
additional smoothness prior. After demonstrating the 
workings of the spectral smoothness prior, we use the 
Gibbs free energy formalism to derive filter formulas for 
the log-normal case and transfer the spectral smoothness 
results to this case in Sec. |III| We demonstrate the per- 
formance of our log-normal reconstruction algorithm in a 



variety of test cases in Sec. IIIB and conclude in Sec. IV 



II. 



RECONSTRUCTING GAUSSIAN FIELDS 
WITH SPECTRAL SMOOTHNESS 



First, we need to introduce some basic assumptions 
and notation. We mainly follow the notation that is used 
throughout the literature on information field theory [e.g. 

El. 

Throughout the paper, we assume that we are analyz- 
ing a set of data d — (c?i, . . . ,dr) € R"" that depends on an 
underlying signal field s : 7W — > R, subject to additive 
noise n = (ni, . . . ,nr) G R"^, 



(1) 



be linear, described by a response operator R, so that 



d = Rs 



(2) 



In most applications, the response operator will include 
some instrumental or observational effects such as obser- 
vations in specific locations of A4 , convolutions of s with 
an instrumental response function, or a Fourier transform 
of the signal field. The only restriction that we make 
here is that the operation that generates the data from 
the signal has to be linear in the signal. Finally, we re- 
strict ourselves to real- valued signals only in the interest 
of notational simplicity. All our results can be straight- 
forwardly generalized to complex-valued signal fields. 
For the noise n, we assume Gaussian statistics, 



(3) 



described by a not necessarily diagonal covariance ma- 
trix N. Here, G{(f>, ^) denotes a multi-variate Gaussian 
probability distribution function with covariance 4>, 



1 



|27r$ 
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cxp 



(4) 



We use the f-symbol to denote a transposed (and in gen- 
eral complex conjugated) quantity and take the limit of 
an infinite-dimensional Gaussian whenever the argument 
is a continuous field, with the appropriate limit of the 
scalar product 



(j)^tlj= / dx (j){x) ipix) y 4>,ip : M 



(5) 



M 



In this section we will deal with the case that the signal 
field s can also be regarded - or at least approximated 
- as a zero-mean Gaussian random field with covariance 
S, 



(6) 



It can be straightforwardly shown [e.g. [7] that the poste- 
rior mean m of the signal field under these assumptions 
is given by 



m = jvssV{s\d) = Dj. 



(7) 



Here, J Vs denotes an integral over the configuration 
space of all possible signal realizations. The operator 



(8) 



is known as the information propagator and the field 

j ^ R^N-^d (9) 



Here, M. is the discrete or continuous space or manifold 
that the signal is defined on. In Sec. |III| we will discuss 
the examples M. = , , and T^. In this section, we 
assume the relationship between signal field and data to 



is called information source. 

In these formulas, the presence of the Gaussian sig- 
nal prior described by the signal covariance serves as a 
means of regularization of the sought continuous signal 
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field reconstruction which otherwise is under-determined 
when only constrained by the finite and noisy data set. 
However, in most physical applications, the signal covari- 
ance is not known a priori. The problem of reconstructing 
Gaussian fields with unknown signal covariances has in 
principle been solved O IHl [131 123 , and even an unknown 
noise covariance can be overcome f21'. 

Enfilin & Weig [S] use the formalism of minimum Gibbs 
free energy to derive filter formulas to be applied to the 
data set when the signal's covariance is unknown. We 
will recapture this formalism briefly in Sec. |III[ where we 
employ it to reconstruct log-normal signal fields. 

Under the assumption of statistical homogeneity and 
isotropy, the unknown signal covariance becomes diago- 
nal in the harmonic eigenbasis, i.e. the Fourier basis for 
signals defined on Euclidean space and the spherical har- 
monics basis for signals defined on the sphere. In the 
following, we denote as k the vector of parameters de- 
termining one mode in the harmonic decomposition, i.e. 
k = (fci, . . . , fc„) for n-dimensional Euclidean space and 
k ~ (i, m) for the two-sphere, where £ is the angular 
momentum quantum number and m the azimuthal one. 
Furthermore, k shall st and for the sc ale of the harmonic 
component, i.e. k — y/kf + ■ ■ ■ + k^ and fc = ^ for R" 
and 5^, respectively. 

Due to the isotropy assumption, the diagonal of the 
signal covariance matrix, the signal's power spectrum P^, 
depends only on the scale k and not on the specific mode 
fc, i.e. 



(10) 



Under these symmetry assumptions, the filter formulas 
derived by EnlBlin & Weig [8] can be written as 



Pk = 



OLk 



Pk 
2 



qk 



2 ^ (i"'^ 

{A;'|fe' = fe} 



mt,l 



Dpk' 



(11) 



(12) 



Note that the first of these equations is the same as in the 
case of a known signal covariance, Eq. ([t]), a generalized 
Wiener filter [7]. Here, ak and qk are parameters used to 
determine the priors for the spectral coefficients Pk (see 
Sec. HI for details) and 



Pk= ^ I 

{k'\k' = k} 



(13) 



is the number of degrees of freedom on the scale k. 

This formalism has been successfully applied for as- 
trophysical reconstructions [THl [20] and in a more ab- 
stract computational setting [24]. In these applications, 
Eqs. (11) and (12 1 were simply iterated. However, they 



were supplemented with an additional ad-hoc smoothing 
step for the power spectrum as part of each iteration. 
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FIG. 1. Power spectra of the one-dimensional test-case. The 
black sohd line shows the theoretical power spectrum, the blue 
dashed line the power in the randomly drawn signal reahza- 
tion studied here, and the green dotted line shows the power 
spectrum reconstructed without smoothing. The horizontal 
dotted line indicates the noise power, given by u^/lOO. 



To illustrate the usefulness of a smoothing step, we 
set up a simple one-dimensional test-case. Here, we as- 
sume that our signal field is defined on the one-sphere 5^, 
i.e. the interval [0, 1) with periodic boundary conditions, 
which we discretize into 100 pixels. Wc make up a power 
spectrum for the signal of the form 



Pk = P 




-7/2 



(14) 



where we choose Pq = 0.2, ko = 5, and 7 = 4. Figures [T] 
and [2] show the power spectrum and the signal realization 
randomly drawn from it, respectively. 

Furthermore, we assume that we have measured the 
signal field in each of the 100 pixels once, subject to ho- 
mogeneous and uncorrelated Gaussian noise with vari- 
ance cr^ = 0.1. In the formalism of Eq. ([2]), this corre- 
sponds to i? = 1 and N = cr^jl. The data realization is 
also shown in Fig. [2] 

We overplot in Figs, [l] and [2] the power spectrum 
and map reconstruction, respectively, obtained by iter- 
ating Eqs. (Ill and ([12]) without smoothing. While the 



signal reconstruction suffers from some too pronounced 
small-scale fluctuations, the reconstructed power spec- 
trum fluctuates wildly for fc > 10 and is therefore not 
trustworthy at all on these scales. The scales on which 
the reconstructed power by chance peaks above the noise 
level are the ones that are especially misrepresented in 
the reconstructed signal field. 

This example should serve to illustrate that some kind 
of spectral smoothing is necessary, especially in cases 
where one is interested in the power spectrum itself. In 
Fig. [2] we also show a comparison to the Wiener filter 
reconstruction, i.e. Eq. ([t]), using the true power spec- 
trum. In the bottom panel of Fig. [2] it can be seen that 
the residuals, i.e. the differences between the true signal 
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FIG. 2. Signal reconstruction in the one-dimensional test- 
case. The blue solid line shows the randomly drawn signal 
realization, the crosses are noisy data points used to do the 
reconstruction, the green dashed line shows the signal recon- 
struction obtained without spectral smoothing and the red 
dotted line is the Wiener filter reconstruction. The lower 
panel shows the same with the true signal substracted. 



A. Reconstruction as a combination of posterior 
mean and maximum a posteriori 

Before incorporating a spectral smoothness prior into 



the derivation of Eqs. (Ill and ( 12 ), we will rederive them 



in a way different from the one presented in [6l [8] . 

As was already mentioned, Eq. ( 11 ) corresponds to the 
posterior mean of the signal under the assumption of a 
known covariance matrix, i.e. a known power spectrum. 
Since the posterior probability distribution is Gaussian in 
this case, its mean also maximizes the probability, so that 
m is the posterior mean and the maximum a posteriori 
solution at the same time, 

m = JvssV{s\d,P) = argmax{V{s\d,P)} . (15) 



Here, we will show that Eq. (12 1 can be written as a 



maximum a posteriori solution as well, considering the 
signal-marginalized posterior 



riP\d) = jvsV{s,P\d). 



(16) 



In order to do this, we first have to define a prior for 
the power spectrum P. In accordance with [BJ [S], we 
choose independent inverse-gamma distributions for each 
spectral component P^, 



and its reconstruction, are reduced if a power spectrum 
is used that is closer to the true one. Smoothness is one 
aspect of the true power spectrum that can be used to 
constrain its reconstruction. 



The two-point correlation function of the signal pro- 
vides another way of looking at the smoothness property 
of its power spectrum. For a statistically homogeneous 
signal, the power spectrum is the Fourier transform of 
the correlation function and vice versa. Thus, a power 
spectrum that exhibits fluctuations on arbitrarily small 
scales in fc-space corresponds to a signal that exhibits 
correlations over arbitrarily large scales in position space. 
Turning this argument around, any signal with correla- 
tions only over a finite range in position space or at least 
correlations that are decreasing rapidly with distance will 
have a power spectrum that does not exhibit any features 
on arbitrarily small scales in fc-space. Thus, the power 
spectrum can be expected to be smooth on the reciprocal 
scale of the typical correlation length scale. 

In the remainder of this section, we will show how to 
incorporate a prior enforcing spectral smoothness into 
the formalism presented thus far. 



n 



1 fPk 

qk^{oLk - 1) \ qk 



exp 



qk_ 



(17) 



Here, r(-) denotes the gamma function, qk is a parame- 
ter defining the location of an exponential cut-off at low 
values of Pk, and ak defines the slope of a power-law for 
large values of Pk- By tuning these parameters, the prior 
can be narrowed or widened according to the a priori 
knowledge about the power spectrum. Taking the limit 
Qk — ^ and Qffc — > 1 turns the inverse-gamma distribution 
into Jeffreys prior which is flat on a logarithmic scale. In 
the examples presented in this paper we always take this 
limit in the final filter formulas. 

In the following, we will work with the logarithmic 
spectral components 



Pk = log Pk 



(18) 



The corresponding prior for these can be straightfor- 
wardly derived from the conservation of probability un- 
der variable transformations and reads 



V{p) = V{P) 

=n 



dP 



dp 



r(afc - 1) 



-(afc-l)pfc 



exp {-que 



Pk ^ 



(19) 
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Using this prior, we can calculate the signal-marginalized 
probability of data and power spectrum, V{d,p), and its 
negative logarithm, the Hamiltonian 

Hid,p) = -logVid,p) 

= - log Q[d - Rs, N) Qis, S) V{p) 

= ^tT {logS)-^tvilogD)-^j^Dj 



where we have made use of the definitions (|8| and ^ and 
have collected all terms that do not depend on the power 
spectrum into an unimportant additive constant. Using 
the spectral dependence of the signal covariance ma- 
trix in a statistically homogeneous and isotropic setting, 
Eq. (10), we can take the derivative of the Hamiltonian 



with respect to one pk and equate it to zero, thus maxi- 
mizing the posterior probability of the logarithmic power 
spectrum. The resulting equation is exactly Eq. (12 1 for 
Pfe e*'* if one makes the identification m 



Dj. 



Thus we have shown that the filter formulas, Eqs. ( 11 ) 



and (12), can be derived as a combination of posterior 



mean for the signal reconstruction and maximum a pos- 
teriori for the power spectrum. This effectively means 
that we have made the approximation 

m = \'Dss'P{s\d) 



(21) 



Vs jVpsV{s\d,p)V{p\d) 

Vs jvpsV[s\d,p)5(p~p^^^^^ 



i.e. we have approximated the posterior probability dis- 
tribution for the power spectrum with a delta distribu- 
tion centered on its maximum. 

It may be worth noting that in the formalism of the 
maximum a posteriori solution for the power spectrum, 
it is straightforward to derive a rough uncertainty esti- 
mate as well. The Hessian of the Hamiltonian gives the 
curvature of the posterior probability distribution and its 
inverse can thus be regarded as an uncertainty matrix. 
For the Hamiltonian given in Eq. ( 20 ) we obtain 



d^H{d,p) 



dpkdpk' 



p=p(MAP) 



= at, - 1 



Pk 



'kk' 



{g|g = fc} 
{q'lq' = k'} 



p=p(MAP) 

(22) 



where 5R(-) denotes the real part of a complex number. 



B. Spectral smoothness priors 

Here, we show how to incorporate a spectral smooth- 
ness prior into the formalism developed in the previous 
section. We do this by augmenting the inverse-gamma 
distributions previously assumed as the spectral prior 
with a probability distribution that enforces smoothness 
of the power spectrum, so that 



Vip)=rsrr.ip)l[VlGiPk). 



(23) 



As an example, we choose a smoothness-enforcing prior 
of the shape 



Vsmip) oc exp 



1 



log Pk 



2^2 /d(logA:) , ,,2 



d (log k) 



(24) 



This prior is maximized by any power-law power spec- 
trum and punishes deviations from such a shape. The 
strength of the punishment is determined by the param- 
eter (Tp. Other useful shapes for a smoothness prior could 
contain the first logarithmic derivative, punishing steep 
spectra, or simple derivatives with respect to fc, punishing 
abrupt changes in the power spectrum. To illustrate the 
meaning of such smoothness priors and point out possible 
caveats, we discuss a few specific cases in App. |A] 

The spectral smoothness prior, Eq. ( 24 ) , can be written 
as a Gaussian in p 



logP, 



Vsmip) exp --p^Tp 



1 



(25) 



where the linear operator T includes both the second 
derivative and the scaling constant a^. The detailed form 
of the operator T that we use in our calculations can be 
found in App. [B] 

Introducing this prior corresponds to adding the term 
■^p^Tp to the Hamiltonian in Eq. (20). Taking its deriva- 



tive with respect to one pk and equating the result with 
zero, we now get 



Qk 



2 S{fc'|/c'=fc} 



D 



1 



Pk 

2 



iTp)k 



(26) 



The only point in which this equation differs from 
Eq. (12) is the extra term Tp in the denominator. 



As before, we can calculate the inverse Hessian of the 
Hamiltonian as an approximate uncertainty matrix. The 
Hessian is 



d^H{d,p) 
dpkdpk' 



p—p\- 



1 



-Pk-Pk' 



{g|5 = fe} 
{q'lq'=k'} 



p=p(MAP) 

(27) 
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FIG. 3. Power spectra of the one-dimensional test-case with 
smoothness prior. The black solid line shows the theoretical 
power spectrum, the blue dashed line the power spectrum re- 
constructed with a smoothness prior with = 1000, and the 
green dotted line the reconstruction with a stifFer smooth- 
ness prior with ap = 10. The hatched regions around the 
dashed and dotted lines are the corresponding one-sigma un- 
certainty regions estimated from the inverse Hessian of the 
Hamiltonian. The horizontal dotted line indicates the noise 
level. 



C. Test-cases 

Using the same one-dimensional test-case as shown in 
Figs, [l] and [2j we perform a reconstruction using the 
spectral smoothness prior given in Eq. (24). The re- 



sulting power spectra using a moderate strength for the 
prior with — 1000 and a strict smoothness prior with 
CTp = 10 are shown in Fig. [s] Clearly, the new recon- 
structions are a much better approximation to the shape 
of the theoretical power spectrum than the one shown 
in Fig. [l] Also, as expected, the power spectrum ob- 
tained when using the strict smoothness prior turns out 
smoother and more closely resembles a power law. 

Also shown in Fig. [3] is an uncertainty estimate for the 
reconstructed power spectra, obtained from the Hessian 
of the Hamiltonian. Taking a one-sigma uncertainty es- 
timate for the logarithmic power spectrum components. 



Hid,p) 




1/2 



(28) 



we plot the uncertainty interval for the power spectra 

as Pk — exp ^p^.'^^^'' ± Ap^ . The uncertainty interval of 

the power spectrum estimates calculated in this way can, 
however, be misleading. It should be noted that the un- 
certainty of the power on the different fc-modes is corre- 
lated. To illustrate this, we plot the complete uncertainty 
matrix, i.e. the inverse Hessian of the Hamiltonian, for 
the case with ap = 1000 in Fig. |4] It can be seen from this 
figure that the correlation of the uncertainty on different 
fc-modes is especially large for small scales. Furthermore, 
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FIG. 4. Inverse Hessian of the Hamiltonian with a smoothness 
prior with = 1000. Plotted is the full matrix given by the 
inverse of Eq. ( |27[ ), evaluated at p = p(MAP) rpj^g diagonal 
of this matrix can be translated into the uncertainty interval 
plotted in Fig.jS] however, that does not take into account the 
correlations across different fc-modes that are visible in this 
plot. 



the inverse Hessian is only a rough approximation to the 
uncertainties in the power spectrum estimation, which 
typically underestimates these. 

In all following examples we will use a spectral smooth- 
ness prior with an intermediate stiffness, given by tip = 
100. 

To study the improvement that the spectral smooth- 
ness prior brings for the reconstruction of the signal field, 
we draw 100 different signal realizations from the same 
power spectrum and add 100 different noise realizations. 
For each of these data sets, we perform the full recon- 
struction and then calculate the power of the difference 
between the reconstructed field m and the true signal 
field s. In Fig.[5]we plot this power, averaged over the 100 
realizations, for different reconstruction schemes. Under 
the assumption that the correct power spectrum is known 
a priori, the residual power is essentially given by the 
noise level on scales for which the signal is dominating 
and by the signal power on scales for which the noise is 
dominating. As can be seen from the plot, reconstruct- 
ing the power spectrum without any spectral smoothing 
leads to a significantly increased residual power on the 
noise-dominated scales, while the quality of the recon- 
structions including the smoothness prior is close to the 
one of the Wiener filter reconstructions, showing only a 
slight excess in residual power on the smallest scales. We 
also plot the average residual power for reconstructions in 
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FIG. 5. Residual power between true signals and their recon- 
structions averaged over 100 different realizations. The black 
solid Une shows the theoretical power spectrum ol the signals, 
given by Eq. ([l4| with Po = 0.2, fco = 5, and 7 = 4. The dot- 
ted horizontal line indicates the noise level, corresponding to 
cr^ = 0.1. The remaining lines show the residual power for 
different reconstruction schemes. From top to bottom, these 
are a reconstruction without any spectral smoothing (green 
dotted line), with an ad-hoc convolution of the power spec- 
trum (magenta dash-dotted line), with the smoothness prior 
given by Eq. ( |24[ ) with cTp = 100 (red short-dashed line), and 
using the correct power spectrum (blue long-dashed line). 




which an ad-hoc smoothing step is apphed to the power 
spectrum after each iteration of Eq. ( [T2| ). This is imple- 
mented as a simple convolution with a Blackman window 
of width Ak = 9. As can be seen in Fig. [5j this ad-hoc 
smoothing partly alleviates the problems of the power 
spectrum reconstruction but is clearly outperformed by 
the rigorous application of a smoothness prior. 

Finally, we consider signal fields of the types discussed 
in App. I A 2] and |A 3] i.e. signals with Gaussian and trian- 
gular correlation functions given by Eqs. (A4) and (A6|, 
respectively. As discussed in the appendix, these cor- 
relation functions lead to theoretical power spectra that 
are strongly suppressed by our spectral smoothness prior. 
Here, we investigate how serious this unwanted effect is 
in practice. 

Shown in Figs. [6] and [7] are the results of the reconstruc- 
tion for the power spectra and signal fields, respectively. 
We again use signals defined on an interval of length one, 
which we divide into 100 pixels, and choose a = 0.2 and 
Co = l/(4V27r) for the Gaussian correlation function, 
Eq. ( [Ail ), s-iid L = 0.2 and Co = 0.25 for the triangular 
correlation function, Eq. (A6|. In both cases, the noise 
variance is = 0.1. 

As expected, the features in the power spectra that are 
discouraged by the spectral smoothness prior are not well 
reconstructed. From Fig. |6] it is obvious that in the case 
with a Gaussian correlation function, the reconstructed 
and true power spectra deviate quite strongly on small 
scales. While the true power spectrum drops off rapidly, 
the reconstructed one stays comparatively flat. In the 



FIG. 6. Pow er sp ectra of signals with corr elation functions 
given by Eq. ( A4[ ) (top panel) and Eq. (A6l (bottom panel). 
The black solid lines show the theoretical power spectra, 
the blue dashed lines the power in the random field real- 
ization drawn from the theoretical power spectra, and the 
green dotted lines the reconstructed power spectra using the 
spectral smoothness prior given in Eq. (24 1 with = 100. 
The parameters substituted into Eq. (A4l are a — 0.2 and 
Co = 1/ (4\/27r) and the ones used in Eq. ( aE I are L = 0.2 and 
Co = 0.25. The noise variance in both cases is cr^ = 0.1. The 
hatched area indicates the uncertainty of the reconstructed 
power spectrum estimated from the inverse Hessian and the 
horizontal dotted line the noise level. 



case of a triangular correlation function, the same can be 
said about the finite fc-value where the true power spec- 
trum drops to zero. While the drop is indeed represented 
by the signal realization, as can be seen from the dashed 
line in the lower panel, the reconstructed power spectrum 
stays relatively level. 

However, the features in the power spectra that the fil- 
ter fails to reconstruct are below the noise level, whereas 
the rise of the power spectrum in the bottom panel of 
Fig. [6] on the smallest scales, which is above the noise 
level, is indeed represented in the reconstructed power 
spectrum as well. Thus, the signal reconstruction does 
not suffer significantly, as can be seen in Fig. [7] An accu- 
rate reconstruction of the power spectrum is important 
mainly at the scales on which the signal-response and 
noise are of comparable magnitude. If the signal power is 
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FIG. 7. Signal fields corresponding to the power spectra 
shown in Fig.|6] The top panel shows a signal with a Gaussian 
two-point correlation function, Eq. (A4l, the bottom panel 
shows one with a triangular correlation function, Eq. ( |A6[ ). 
Shown in each plot are the signal realization s (blue solid 
line), the data d (crosses), the signal reconstruction m (green 
dashed line), and an estimate of the local one-sigma uncer- 
tainty of the reconstruction, given by m±diag(D)^''^ (hatched 
area). At the bottom of each panel, the same is plotted with 
the true signal subtracted. 



dominant, the reconstruction will follow the data closely, 
while it will smooth the data heavily in the opposite case 
of dominating noise, irrespective of the exact shape of 
the power spectrum. Thus, the reconstruction algorithm 



will in general perform well even in cases in which cer- 
tain features in the power spectra are not allowed by the 
spectral smoothness prior. However, if the objective is an 
accurate reconstruction of the power spectrum and any 
such features are suspected to be present, the spectral 
smoothness prior needs to be adapted to this situation. 



III. RECONSTRUCTING LOG-NORMAL 
FIELDS 

Now we turn to the problem of reconstructing a log- 
normal field. We define our signal field s to be the loga- 
rithm of this log-normal field, so that the prior probabil- 
ity distribution for s is again a Gaussian, described by a 
mean and a covariance S. For simplicity, we assume the 
prior mean to be zero. The data model then becomes 



d = Rc' 



(29) 



where we have again included additive Gaussian noise n 
and the application of the exponential function to the 
signal field is to be interpreted pointwise. Due to the 
non-linearity of the exponential function, the posterior 
of s. 



r{s\d, s)(xg{d-Re',N)g{s,s), 



(30) 



is highly non-Gaussian, even when the signal covariance 
S is known. Adding a prior for S and marginalizing over 
it only makes the problem more complex, 

V{s\d) (xGid- Re'',N) jvsg{s,S)V{S). (31) 

The path we pursue here is to treat this probability dis- 
tribution approximatively. In order to do so, we employ 
the formalism of minimum Gibbs free energy presented 
in [8]. The basic idea is to approximate the posterior 
V{s\d) with a Gaussian, described by a mean m and a 
covariance D. An approximate Gibbs free energy can be 
calculated as a function of these quantities and it was 
shown in [8] that the optimal Gaussian approximation - 
in the sense of minimum KuUback-Leibler divergence - 
can be obtained by minimizing the approximate Gibbs 
free energy with respect to m and D. 

The approximate Gibbs free energy is defined as 

G{m,D) = tj{m,D)-TSB{D). (32) 
The last term in this definition is the Boltzmann entropy 



^B(^)-^tr(l + log(27ri?)), 



(33) 



which depends only on the covariance D. The other term 
in the approximate Gibbs free energy is the approximate 
internal energy 



U{m,D) = {H{s,d))^^^_^^j,y 



(34) 
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Here, 



ig{s 



-m D) denotes an expectation value calculated 
with respect to the Gaussian posterior approximation 
and H{s,d) — —log'P{s,d) is the Hamiltonian of the 
problem. Finally, the temperature T can be regarded as 
a tuning parameter that regulates the importance that 
is given to the region around the maximum of the full 
posterior and regions removed from the maximum. The 
limiting value T — leads to the maximum a posteriori 
solution for the signal field. In the following, we choose 
the default value of T = 1 and refer the reader to [8lf9l[TT| 
for an in-depth discussion of the temperature parameter. 

Assuming independent inverse-gamma distributions as 
priors for the power spectrum components Pk, as we 
had already done in Eq. (17), the calculation of the 
Hamiltonian yields 

H{s,d) - ~ log {V{d\s)V{s\P)r{P)) 

= - log {g{d-Re',N)gis,s)l[riG{Pk)] 



{k'\k'=k} 



+ const. 



(35) 



Here, we have again collected s-independent terms in an 
additive constant and introduced the abbreviation 



(36) 



In order to calculate the Gaussian expectation value of 
this Hamiltonian analytically, we expand the logarithm 
appearing in this expression in a power series around its 
expectation value 




(37) 



y {k'\k'=k} ) 

log [qk) 




qk 



(38) 



We truncate this expansion after the first order, i.e. 
*max = 1- Note that our choice of q}~ ensures that the 
first-order term itself vanishes. 



With this simplification, we can calculate the approx- 
imate Gibbs free energy to be 



G{m,D) = / dxj^e 
Im 



+ [ dx [ dy-A/4j,e™="+'"«+5^"=- + 5^»»+^== 
Jm Jm 2 



-t- 



'M 

E 

k 



M 

ak-l 



Pk 



log (qk) 



1 



tr(H- log (2771?)). 



(39) 



To avoid confusion, we write out explicitly all integrals 
appearing here and in the following. Taking the func- 
tional derivatives with respect to m and D and equat- 
ing them with zero yields two filter equations that de- 
termine m and D. In these equations, the right hand 
side expression of Eq. (12 1 appears. If we reidentify this 



expression with the spectral components Pk and write 
^kk' ~ ^kk'^f^' filter equations become 



dyS^ 



M 



xy 



jy<d^^^ 2 ^yy 



- / dzMj,^e'"«+'"- + 3^'"'+5^"+-°"- 
Jm 



(40) 



and 



D- 



+ / dzM^^e^^+^^ + s^^-^ + s-D^.+i^..; 



M 



Jxy 



V / xy 



(41) 



Together with Eq. (12 1, the last two equations fully 



determine the Gaussian approximation to the posterior. 
Solving them self-consistently gives an estimate m for 
the signal field and an estimate D for the correspond- 
ing uncertainty matrix. The corresponding approximate 
posterior mean estimate of the exponentiated field values 
then is 



/V{s\d) 



ig{s-m,D) 



(42) 



Of special interest for many applications is the case in 
which the matrix M = R^N~^R is diagonal, e.g. when 
the noise contributions to the individual data points are 
uncorrelated and the response is purely local. In this case 
the filter equations simplify somewhat to 



= dy S^. 



M 



Jy 



(43) 



and 



V / xy 



(44) 
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A. Spectral smoothness in the log-normal case 



In Sec. Ill Al we had seen that in the Gaussian case the 
full reconstruction with unknown power spectrum can be 
regarded as a combination of the posterior mean recon- 
struction under the assumption of a known power spec- 
trum and the estimation of this power spectrum as the 
one that maximizes its posterior probability. In the case 
of a log-normal field, the posterior mean cannot be calcu- 
lated analytically, even if the power spectrum is assumed 
to be known. However, it is interesting to note that em- 
ploying the formalism of minimum Gibbs free energy to 
the log-normal reconstruction problem with known power 
spectrum, one arrives exactly at the formulas derived in 



the previous section, Eqs. (40 1 and (41) 



Thus, the full reconstruction, consisting of Eqs. (401, 



(41 1, and ( 12 1 can again be regarded as a combination of 



the calculation of the posterior mean for the signal under 
the assumption of a power spectrum and the estimation 
of the power spectrum according to Eq. (12 1. From this 



viewpoint, the inclusion of a smoothness prior for the 
power spectrum is trivial. We simply replace the power 
spectrum estimation step according to Eq. ( 12 ) with the 



one derived in Sec. II B i.e. with Eq. (26 1, just as we had 



done in the Gaussian case. 




80 100 



B. 



Test cases 



In this section, we present some test cases for the the- 
ory developed so far on the reconstruction of log-normal 
fields. We study one-dimensional tests with differing de- 
grees of non-linearity and differing noise-levels, as well 
as two two-dimensional test-cases. For simplicity, we as- 
sume an ideal local response, i? = 1, except in the last 
example, where we study the effects of an observational 
mask. The noise is assumed to be uncorrelated and ho- 
mogeneous, N — cr^l. F or the smoothness prior, we use 

and choose — 100. 



II B 



the one discussed in Sec. 

First, we discuss again tlie case of a signal on the one- 
sphere S^, i.e. the interval [0, 1) with periodic boundary 
conditions, here discretized into 100 pixels. Shown in 
Figs. [8p0| are three different cases. In each of these cases, 
the signal realization is the same as in Sec.|llj only with a 
differing normalization of the power spectrum. In Fig. [8] 
we show a mildly non- linear case with Pq = 0.1 and in 
Figs. |9] and [To] a highly non-linear case with Pq = 0.3. 
The noise variance cr^ is 0.1 in Fig. [s] 1 in Fig. [9] and 25 
in Fig. [To] In the lower panels of these three figures, we 
show the signals and their reconstructions and in the up- 
per panels the exponentiated fields, e*, as well as our pos- 
terior mean estimate for these, (e'')p(g|^j « g"i-i-2diag(r>)_ 
In the mildly non-linear case the reconstruction is a 
reasonably good approximation to the true signal both 
in regions of high signal values and regions of low signal 
values. However, it is apparent that the quality of the 
reconstruction is slightly higher in the former regions. 
This is due to the homogeneity of the noise statistics that 



FIG. 8. One-dimensional log-normal reconstruction in a 
mildly non-linear low-noise case. The power spectrum from 



which the signal realization was drawn is given by Eq. 1 14 1 
with Po — 0.1, fco ~ 5, and 7 = 4. The noise variance is 
= 0.1. The top panel shows the exponentiated signal field 
e" (blue solid line), the data d (crosses), the reconstruction 
gm+jdiag(D) ('gj.ggjj dashed line), and the uncertainty inter- 
val of the reconstruction, given by e'"+5<i'='s(D)±{diag(D))i/2 
(hatched region). The lower panel shows the signal field s 
(blue solid line), the logarithm of the data log(ci) (crosses), 
the reconstruction m (green dashed line), and the uncer- 
tainty interval for m, given by m± (diag(Z)))^''^. In the lower 
panel, only data points for which logd is greater than —1.5 
are shown. 



we have assumed. Since our signal response, given by 
-Re", depends non- linearly on the signal, the noise impact 
is lower in regions of higher signal values and hence the 
signal inference is less demanding in these regions. From 
Fig. [8] it is apparent that this effect is well represented by 
the point-wise uncertainty of the reconstruction, given by 
diag(D). As can be seen from Fig. [o] the effect becomes 
more pronounced for higher degrees of non-linearity, i.e. 
larger signal variances on a logarithmic scale. 

The highly nonlinear case with high noise level, de- 
picted in Fig. 
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exhibits a point-wise uncertainty esti- 
mate for the reconstruction that can clearly not be in- 
terpreted as a 68% confidence interval. This is due to a 
known problem of the filter discussed in Sec. JH] and by 
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FIG. 9. One-dimensional log-normal reconstruction in a 
highly non-linear low-noise case. The plotted quantities are 
the same as in Fig. |8] The parameters describing the power 
spectrum are Po = 0.3, fco = 5, and 7 = 4 and the noise vari- 
ance is £7^ = 1. In the lower panel, only data points for which 
logd is greater than —3 are shown. 
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FIG. 10. One-dimensional log-normal reconstruction in a 
highly non-linear high-noise case. The plotted quantities are 
the same as in Fig. [S] The parameters describing the power 
spectrum are Po ~ 0.3, fco = 5, and 7 = 4 and the noise 
variance is = 25. In the lower panel, only data points for 
which logd is greater than —3 are shown. 



extension also of the filter derived in this section. As dis- 
cussed in [6], the filter exhibits a perception threshold. 
This means that if the signal-response-to-noise ratio is 
lower than a certain threshold on a given scale, then the 
filter will not reconstruct any power on this scale. Our 
usage of the spectral smoothness prior partly alleviates 
this porblem in that it prevents the power of individ- 
ual scales to drop to zero. However, the reconstructed 
power on the noise-dominated scales will in general still 
be too low. This directly affects the estimate for the re- 
construction's uncertainty, given hy D = (^S~^ -f Af) , 
which tends toward zero in the limit of zero power on all 
scales. 

Furthermore, the lack of power on all but the few 
signal-dominated scales can lead to ringing effects, i.e. 
prominent signal-dominated features are well recon- 
structed and extrapolated periodically into the noise- 
dominated regions. The deep trough in the signal recon- 
struction that can be seen in the lower panel of Fig. 10 



around pixel number 35 is most likely due to this effect. 
The high degree of non-linearity acts to reinforce this 
effect. As can be seen in the top panel of Fig. 10 the ef- 



fect of the trough onto the exponentiated reconstruction 
is almost intangible. 

It is certainly mandatory to keep the potential prob- 
lems of ringing and an underestimated uncertainty in 
mind when reconstructing a field which is swamped by 
noise in the better part of its domain. However, as can 
be seen from Fig. [Toj the main features of the signal field 
will still be reconstructed reliably even in such an unfa- 
vorable case. 

In these test cases, we have chosen the one-dimensional 
interval as domain of the log-normal signal field mainly 
for illustrative purposes. The algorithm is, however, ver- 
satile in that it can operate on virtually any spacaM To 
illustrate this point, we close this section with two two- 
dimensional examples. 

In Fig. [TT] wc show an example for a signal defined 
on the two-torus = x S^, i.e. a segment of 



^ We are using the NIFTY package |23| in our implementation. This 
makes changing the domain of the signal field trivial, requiring 
only minuscule changes to the code. 
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FIG. 11. Log-normal reconstruction example in a toroidal setting. The signal's power spectrum is given by Eq. (141 with 
Po ~ 0.3, fco — 2, and 7 = 5. The noise variance is = 10. The top row shows the data set d (left) and its logarithm (right). 
In the logarithmic version, pixels with negative data values are plotted as dark blue. The middle row shows the exponentiated 
signal field e" (left) and its reconstruction, given by e'""'"2'*"'^'^^\ (right). The corresponding non-exponentiated quantities s 
(left) and m (right) are shown in the bottom row. 
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FIG . 12 . Log-normal reconstruction example on a two-sphere. The quantities plotted in the top three rows are the same as in 
Fig. |ll| The bottom row shows the uncertainty of the signal reconstruction given by diag(_D)^''^ (left side) and the fractional 

uncertainty of the reconstruction of the exponentiated field given by e'*'^^'^' ^ —1 (right side). The power spectrum parameters 
are Po ~ 0.3, kg — 5, and 7 = 4 and the noise variance is ct^ — 10. In the top row, pixels without measurements are plotted 
gray. 
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with periodic boundary conditions. The periodicity here 
and in the earlier one-dimensional examples is a conse- 
quence of our use of Fast Fourier Transform routines. It 
can be approximately overcome by artificially extending 
the interval or the section of for more than the typi- 
cal correlation length of the signal while leaving the data 
unchanged when doing the reconstruction. Thus, the in- 
fluence of the reconstruction on one side of the interval 
or rectangle will have negligible influence on the recon- 
struction on the other side. 

We discretize the two-torus and its corresponding 
Fourier space into 50 x 50 pixels. The rectangular orien- 
tation of the Fourier pixels leads to many different scales 
k being represented, each, however, only by a few pix- 
els, corresponding to a few different Fourier vectors k. 
To avoid dealing with all these scales individually, we 
bin the scales logarithmically into 17 bins. We replace 
all summations over Fourier components with a certain 
scale appearing in the filter formulas with summations 
over all Fourier components whose scale falls within one 
bin. Thus, we reconstruct the power for each bin instead 
of for each scale that is represented in the rectangular 
Fourier grid. The signal is again drawn from the power 
spectrum given by Eq. (14) with Pq = 0.3, ko = 2, and 
7 = 5. The noise level is chosen as = 10. Signal, data, 
and reconstruction are shown in Fig. |11[ 

The toroidal example exhibits essentially the same fea- 
tures that we had already seen in the one-dimensional 
case. While the reconstruction is generally in good agree- 
ment with the true underlying signal, it is less accurate 
in the regions of small signal values where the signal- 
response-to-noise ratio is far worse due to the quite high 
degree of non-linearity. 

In Fig. \12\ we show an example for a reconstruction on 
the two-sphere S^. This example has special relevance 
for astronomical applications since one can interpret the 
sky as a two-sphere and therefore any astrophysical sig- 
nals without distance information will be defined on this 
manifold. We generate a mock signal from the power 
spectrum given in Eq. ( [T4| with Pq — 0.3, ko = 5, and 
7 = 4. The noise level is = 10. We discretize the 
sphere using the HEALPixrl package [10] with a reso- 
lution parameter of A^side = 16 leading to 3 072 pixels 
in total. The power spectrum components Pk are in this 
case the components of an angular power spectrum, often 
denoted as C^. 

In this last case study, we replace the trivial response 
R = 1 with a projection onto part of the sphere, effec- 
tively masking 600 pixels around the equator for which 
we assume that no measurements have been taken. This 
resembles a typical situation in extragalactic astronomy, 
where observations through the Galactic plane are not 
possible due to the obscuration by the Milky Way. 



^ The HEALPlx package is available from |http : //healplx ■ jpl . | 
I" nasa.gov/^ . 



Fig. 12 shows the signal field, data, and reconstruc- 
tion both in the exponentiated and in the linear version. 
In the panels showing the reconstructed signal field, it 
can be nicely seen how the algorithm is able to extrapo- 
late into the gap region from the data on the boundary. 
This is possible due to the knowledge of the correlation 
structure that was inferred from the same data set. Also 
shown in Fig. [T2] are the pixel-wise uncertainty estimate 

1/2 

of the signal field's reconstruction, given by diag(I?) , 
and the fractional uncertainty of the reconstructed expo- 
nentiated signal, given by e'^"*s(£')^''^ _ which is approx- 

1 /2 

imated well by diag(£') in most regions. It can be seen 
that the uncertainty tends to be higher in the regions of 
low signal values, as was the case in the one-dimensional 
examples, and in the region around the equator which 
lacks observations. This is to be expected, since the only 
constraint on the signal in this region comes from ex- 
trapolations from neighboring regions using the signal's 
correlation structure inferred from the data. 



IV. SUMMARY AND CONCLUSIONS 

We have developed an algorithm to infer log-normal 
random fields from noisy measurement data. The log- 
normal model was chosen due to its wide range of applica- 
tions in astrophysics and other fields. The reconstruction 
method uses the correlation structure of the log-normal 
field to differentiate between features in the data that 
are due to noise and such that are due to variations in 
the true underlying field. This correlation structure, de- 
termined by the field's power spectrum, is, however, in 
general not known a priori. We have therefore extended 
the theory for simultaneous reconstructions of a field and 
its power spectrum that was developed and applied suc- 
cessfully for Gaussian random fields in the past [6l |8j [T9l - 
EU [24] to log-normal fields. 

An additional feature of our reconstruction method is 
the use of a smoothness prior for the power spectrum. 
We have suggested to employ a prior based on the sec- 
ond double-logarithmic derivative of the power spectrum 
and shown that it is well suited to handle a large variety 
of cases. Having investigated possible pitfalls associated 
with the usage of such a prior from theoretical as well as 
practical viewpoints, we should stress that the derivation 
of the filter formulas laid out here does not depend on the 
specific form of the spectral smoothness prior. In cases 
in which the prior we employed here cannot be expected 
to yield satisfactory results it should simply be replaced 
by a different one. 

The algorithm we have derived depends in no way on 
the space on which the signal field is to be reconstructed. 
We have demonstrated this by showing examples of re- 
constructions of mock signals on a one-dimensional inter- 
val, a flat two-dimensional space, and a spherical space. 
We have discussed the performance of the algorithm in 
these scenarios and pointed out possible caveats when 
dealing with very low signal-response-to-noise ratios. 
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In these application examples, we have assumed that 
the observational data simply represent the underlying 
log-normal field, subject to additive noise. Furthermore, 
we have illustrated the ability of the algorithm to extrap- 
olate from the given data in a test case in which these 
data were assumed to be incomplete, thus demonstrat- 
ing the power of the usage of the correlation information 
contained in the data. The derivation of the filter formu- 
las, however, is even more general. It allows for an ar- 
bitrary linear relationship between log-normal field and 
data, described by a response matrix R. The resulting 
formulas include the general response matrix which can 
e.g. represent an incomplete observation, a convolution, 
or a Fourier transformation. We have also allowed for a 
general noise covariance matrix, thus including cases of 
heteroscedastic or correlated noise. 



This makes the algorithm widely applicable. 
Applications that we have in mind include for ex- 
ample the study of diffuse Galactic emission components 
at radio frequencies and reconstructions of emissiv- 
ity fields across galaxy clusters from interferometric 
observations. 
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for example for Brownian motion of a particle, and of an 
exponential two-point correlation function, i.e. 



C{t) = S^y = Coe''^^ 



(A3) 



where r — |a: — y| is the distance between two points or 
time-instances x and y. Such a correlation function arises 
for example from the Ornstein- Uhlenbeck process |25j . 
Fourier transforming this correlation function in a one- 
dimensional space yields the power spectrum which takes 
on the form given by Eq. (14 1 with Pq — fco = /3, 
and 7 = 2. 

However, there are also cases that do not lead exclu- 
sively to small values of the double-logarithmic deriva- 
tives. Two such scenarios will be studied in the remain- 
der of this appendix. 



Appendix A: Case studies for the spectral 
smootiiness prior 

In this appendix, we discuss a few possible correlation 
structures and power spectra and examine the suitability 
of a spectral smoothness prior of the form described by 
Eq. pi. 



1. Power law spectra 



If the signal's power spectrum is a broken power law of 



the shape given by Eq. (14), the first double- logarithmic 
derivative is 



dpk 
dlogk 



7 



(Al) 



which behaves like —7 

ward —7 for k — > 00. 
derivative is 

d^Pk 



(^-^^ for fc <C fco and tends to- 
The second double-logarithmic 



27 



d (log k) 



2 ' 



(A2) 



which tends toward zero both for fc — > and fc — )■ 00. The 
second derivative takes on its extremum, given by — ^, 
at the knee frequency fcg where the spectral index of the 
power spectrum changes from to 7. Thus neither a prior 
punishing large values for the first logarithmic derivative, 
nor one punishing large values for the second logarithmic 
derivative prevents the true solution from being found, 
provided the values of —7 and — ^, respectively, are well 
within the range of values that are allowed for the deriva- 
tive by the prior. In fact, choosing Up = 7/2 would allow 
such a change in spectral index roughly once per e-folding 
of the fc-value. Choosing fp = 57^ would turn such a kink 
into a less common /i-sigma event. 

The case of a broken power-law contains the special 
cases of a pure power-law, i.e. fco — ^ 0, such as arises 



2. Gaussian correlations 

Consider a two-point correlation function of Gaussian 
shape for a field on a one-dimensional space, i.e. 



C{r) = Coe- 



(A4) 



A field whose statistics are described by such a two- 
point correlation function can be regarded as a station- 
ary and spatially uncorrelated field convolved with a 
Gaussian. Calculating the corresponding power spec- 
trum via Fourier transformation yields 



(A5) 



This power spectrum drops quickly with increasing fc, due 
to the flatness of the correlation function around r — 0. 
Therefore, both the first and second double-logarithmic 
derivatives grow unbounded as fc — > 00. Thus, by employ- 
ing a spectral smoothness prior that punishes large values 
for these derivatives, one prevents in principle the recon- 
struction of the true power spectrum and supresses the 
small-scale correlations in the reconstructed field that are 
in reality more pronounced. In cases in which Gaussian 
correlations are expected, it might therefore be advisable 
to choose the strength of the spectral smoothness prior, 
given by (Tp, fc-dependent or choose a different smooth- 
ness prior altogether. 



3. Triangular correlations 

Another case in which the double-logarithmic deriva- 
tives of the power spectrum can become divergent is a 
correlation function with finite support. As an example, 
we consider a signal on a one-dimensional space with cor- 
relations only over the finite distance 2L, given by 



Cir 







for r < i 
else 



(A6) 
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This triangular correlation function describes a station- 
ary and spatially uncorrelated field that has been con- 
volved with a top-hat kernel. Fourier transforming it 
yields the power spectrum 



Pfe-^(l-cos(fcL)). 



(A7) 



This power spectrum becomes exactly zero at finite k- 
values. Its logarithm, and therefore also the double- 
logarithmic derivatives, are divergent at these locations, 
so not even a fc-dependent value of tXp can ensure the cor- 
rect reconstruction of the power spectrum in this case. 
The oscillatory behavior of this power spectrum is a 
generic feature of signal fields that are correlated only 
over a finite distance. 

Note that both in the case of only locally correlated 
fields and in the case of Gaussian correlations, the double- 
logarithmic derivatives can be kept finite by adding a 
constant floor to the power spectrum, i.e. introducing 
an additive part to the signal that is spatially uncorre- 
lated. In many cases, the variance of the uncorrelated 
addition needed to satisfy the spectral smoothness prior, 
Eq. ( 24 1 , will be small enough so as not to influence the 



signal reconstruction significantly. Note that in practice, 
any field variations are restricted by the finite pixel size 
and only fc-values up to a finite fcmax will be considered. 
In Sec. |II C[ we investigate the effect that the spectral 
smoothness prior given in Eq. ( 24 1 has in practice on 



the recunstruction of a signal field exhibiting the two po- 
tentially problematic two-point correlations that we dis- 
cussed here. 

In conclusion, using a spectral smoothness prior that 
punishes large values for the first or second double- 
logarithmic derivative of the power spectrum can lead 
to the introduction of spurious small-scale variations in 
cases in which the true two-point correlation function is 
flat around r = or has only finite support. In the latter 
case, it can also introduce spurious large-scale correla- 
tions. This will, however, in general only be a problem 
for the reconstruction if some feature that mimics these 
large-scale correlations is present in the data, i.e. caused 
by the noise. Another case in which the employment of 
any spectral smoothness prior is obviously a bad idea 
is that of a signal that exhibits prominent periodicities. 
The detection of such spectral lines would only get hin- 
dered by the usage of a spectral smoothness prior. 



Appendix B: Discretization of the spectral 
smoothness prior 

In our implementation we use discretized values 
(^i)i=o i„ax length scales that are represented 



on the computational grid, or in case of the two-torus, 
bins thereof. We approximate the integral in the expo- 
nent of the spectral smoothness prior, Eq. ( 24 1 , with a 
sum according to 



d(log.) 



J2 5d^p)l (Bl) 



Here, we use the abbreviations 

log/Ci+i - logfci_i 

0,- = 

2 

for integer values of i and 

5^ = logfcj+i/2 - l0gfcj_i/2 



(B2) 



(B3) 



for half-integer values of i. Note that we have excluded 
the boundaries at and ki_^^^ from the sum to avoid 
numerical problems at these locations. 

We approximate the second logarithmic derivative as 



3 * 



SO that we can represent it as a matrix with the entries 



Ai^iii 



- ( ^ I ^ 

V ^4+1/2 <^i-l/2 

1 



(B5) 
(B6) 



acting on the vector p. All other entries of the matrix A 
are zero. 



We can now write the exponent of Eq. ( 24 1 as 



o 2 Aaogfc) f^^^V^-^E^^^^ 



where the matrix T is given by 



(B7) 
(B8) 



